      Subroutine Aver_Ghost(FIELD, GX, GHOSTDIR)
      include "param.fi"
      include "common.fi"
      REAL FIELD(NTOT3T), GX(NTOT3T)
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      Real GHOSTI, GHOSTJ, GHOSTK, FIB2(3)
      INTEGER GHOSTDIR

      Do IB = 1, NTOT3T
        GX(IB) = FIELD(IB)
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        GHOSTI = 1.0
        I = HASHIND_Ghost(II + DELI)
        If (I .LT. 0) Then
          GHOSTI = -1.0
          I = ABS(I)
        Endif
        Do DELJ = -1, 1
        GHOSTJ = 1.0
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND_Ghost(JJ + DELJ)
        If (J .LT. 0) Then
          GHOSTJ = -1.0
          J = ABS(J)
        Endif
        Do DELK = -1, 1
        GHOSTK = 1.0
        K = HASHIND_Ghost(KK + DELK)
        If (K .LT. 0) Then
          GHOSTK = -1.0
          K = ABS(K)
        Endif
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC_Ghost(X1(IB2),X2(IB2),X3(IB2),
     &                                X1(IB),X2(IB),X3(IB),
     &                                GHOSTI,GHOSTJ,GHOSTK,
     &                                TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX + TMPY*TMPY + TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance.GT.1E-10
     &         .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)
                FIB2(1) = FIELD(IB2)
                If (GHOSTI.LT.0 .AND. GHOSTDIR.EQ.1) FIB2(1) = -FIB2(1) 
                If (GHOSTJ.LT.0 .AND. GHOSTDIR.EQ.2) FIB2(1) = -FIB2(1) 
                If (GHOSTK.LT.0 .AND. GHOSTDIR.EQ.3) FIB2(1) = -FIB2(1) 
              If (Shepard.EQ.0) Then
                GX(IB) = GX(IB) + 0.5* (FIB2(1)-FIELD(IB)) * Weight(IB2)
     &                 * Kernel(Distance)
              Else
                GX(IB) = GX(IB) + 0.5* (FIB2(1)-FIELD(IB)) * Weight(IB2)
     &                 * Kernel(Distance) / SheCoef(IB)
              Endif  
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine
C==================================================================
      Subroutine Grad_Ghost(FIELD, GX, GY, GZ, GHOSTDIR)
      include "param.fi"
      include "common.fi"
      REAL FIELD(NTOT3T), GX(NTOT3T), GY(NTOT3T), GZ(NTOT3T)
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      Real GHOSTI, GHOSTJ, GHOSTK, FIB2(3)
      INTEGER GHOSTDIR

      Do IB = 1, NTOT3T
        GX(IB) = 0.0
        GY(IB) = 0.0
        GZ(IB) = 0.0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        GHOSTI = 1.0
        I = HASHIND_Ghost(II + DELI)
        If (I .LT. 0) Then
          GHOSTI = -1.0
          I = ABS(I)
        Endif
        Do DELJ = -1, 1
        GHOSTJ = 1.0
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND_Ghost(JJ + DELJ)
        If (J .LT. 0) Then
          GHOSTJ = -1.0
          J = ABS(J)
        Endif
        Do DELK = -1, 1
        GHOSTK = 1.0
        K = HASHIND_Ghost(KK + DELK)
        If (K .LT. 0) Then
          GHOSTK = -1.0
          K = ABS(K)
        Endif
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC_Ghost(X1(IB2),X2(IB2),X3(IB2),
     &                                X1(IB),X2(IB),X3(IB),
     &                                GHOSTI,GHOSTJ,GHOSTK,
     &                                TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX + TMPY*TMPY + TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance.GT.1E-10
     &         .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)
                FIB2(1) = FIELD(IB2)
                FIB2(2) = FIELD(IB2)
                FIB2(3) = FIELD(IB2)
                If (GHOSTI.LT.0 .AND. GHOSTDIR.EQ.1) FIB2(1) = -FIB2(1) 
                If (GHOSTJ.LT.0 .AND. GHOSTDIR.EQ.2) FIB2(2) = -FIB2(2) 
                If (GHOSTK.LT.0 .AND. GHOSTDIR.EQ.3) FIB2(3) = -FIB2(3) 
              If (Shepard.EQ.0) Then
                GX(IB) = GX(IB)-KX *(FIB2(1)-FIELD(IB)) * Weight(IB2)
                If (DIMM.EQ.DIMM3D)
     & GY(IB) = GY(IB)-KY *(FIB2(2)-FIELD(IB)) * Weight(IB2)
                GZ(IB) = GZ(IB)-KZ *(FIB2(3)-FIELD(IB)) * Weight(IB2)
              Else
                GX(IB) = GX(IB)-KX *(FIB2(1)-FIELD(IB)) * Weight(IB2)
     &                 / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     & GY(IB) = GY(IB)-KY *(FIB2(2)-FIELD(IB)) * Weight(IB2)
     &                 / SheCoef(IB)
                GZ(IB) = GZ(IB)-KZ *(FIB2(3)-FIELD(IB)) * Weight(IB2)
     &                 / SheCoef(IB)
              Endif  
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine
C==================================================================

      Subroutine Calc_ShepardCoef_Ghost
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      Real GHOSTI, GHOSTJ, GHOSTK, FIB2(3)

      Do IB = 1, NTOT3T
        SheCoef(IB) = 0.0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        GHOSTI = 1.0
        I = HASHIND_Ghost(II + DELI)
        If (I .LT. 0) Then
          GHOSTI = -1.0
          I = ABS(I)
        Endif
        Do DELJ = -1, 1
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        GHOSTJ = 1.0
        J = HASHIND_Ghost(JJ + DELJ)
        If (J .LT. 0) Then
          GHOSTJ = -1.0
          J = ABS(J)
        Endif
        Do DELK = -1, 1
        GHOSTK = 1.0
        K = HASHIND_Ghost(KK + DELK)
        If (K .LT. 0) Then
          GHOSTK = -1.0
          K = ABS(K)
        Endif
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC_Ghost(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          GHOSTI,GHOSTJ,GHOSTK,
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP) Then
                SheCoef(IB) = SheCoef(IB) +
     &                        Kernel(Distance) * Weight(IB2)
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine

C==================================================================

      Subroutine Calc_dRho_dT_Ghost
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      Real GHOSTI, GHOSTJ, GHOSTK, FIB2(3)

      Do IB = 1, NTOT3T
        Temp1(IB) = 0.0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        GHOSTI = 1.0
        I = HASHIND_Ghost(II + DELI)
        If (I .LT. 0) Then
          GHOSTI = -1.0
          I = ABS(I)
        Endif
        Do DELJ = -1, 1
        GHOSTJ = 1.0
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND_Ghost(JJ + DELJ)
        If (J .LT. 0) Then
          GHOSTJ = -1.0
          J = ABS(J)
        Endif
        Do DELK = -1, 1
        GHOSTK = 1.0
        K = HASHIND_Ghost(KK + DELK)
        If (K .LT. 0) Then
          GHOSTK = -1.0
          K = ABS(K)
        Endif
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC_Ghost(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          GHOSTI,GHOSTJ,GHOSTK,
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP .AND. Distance.GT.1E-10
     &         .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)
                FIB2(1) = U1(IB2)
                FIB2(2) = U2(IB2)
                FIB2(3) = U3(IB2)
                If (GHOSTI.LT.0) FIB2(1) = -FIB2(1) 
                If (GHOSTJ.LT.0) FIB2(2) = -FIB2(2) 
                If (GHOSTK.LT.0) FIB2(3) = -FIB2(3) 
              If (Shepard.EQ.0) Then
                Temp1(IB) = Temp1(IB) -
     &            KX * (FIB2(1)-U1(IB)) * Weight(IB2) * Rho(IB2)
                If (DIMM.EQ.DIMM3D)
     &           Temp1(IB) = Temp1(IB) -
     &            KY * (FIB2(2)-U2(IB)) * Weight(IB2) * Rho(IB2)
                Temp1(IB) = Temp1(IB) -
     &            KZ * (FIB2(3)-U3(IB)) * Weight(IB2) * Rho(IB2)
              Else
                Temp1(IB) = Temp1(IB) -
     &            KX * (FIB2(1)-U1(IB)) * Weight(IB2) * Rho(IB2)
     &            / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     &           Temp1(IB) = Temp1(IB) -
     &            KY * (FIB2(2)-U2(IB)) * Weight(IB2) * Rho(IB2)
     &            / SheCoef(IB)
                Temp1(IB) = Temp1(IB) -
     &            KZ * (FIB2(3)-U3(IB)) * Weight(IB2) * Rho(IB2)
     &            / SheCoef(IB)
              Endif
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine
C==================================================================

      Subroutine Calc_dU_dT_Ghost
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      INTEGER tagGHOST
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      REAL TMPXR, TMPYR, TMPZR, Distance_Real
      Real GHOSTI, GHOSTJ, GHOSTK, PIB2

      Do IB = 1, NTOT3T
        Temp1(IB) = 0.0
        Temp2(IB) = 0.0
        Temp3(IB) = 0.0
        tagGHOST = 0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        GHOSTI = 1.0
        I = HASHIND_Ghost(II + DELI)
        If (I .LT. 0) Then
          GHOSTI = -1.0
          tagGHOST = 1
          I = ABS(I)
        Endif
        Do DELJ = -1, 1
        GHOSTJ = 1.0
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND_Ghost(JJ + DELJ)
        If (J .LT. 0) Then
          GHOSTJ = -1.0
          tagGHOST = 1
          J = ABS(J)
        Endif
        Do DELK = -1, 1
        GHOSTK = 1.0
        K = HASHIND_Ghost(KK + DELK)
        If (K .LT. 0) Then
          GHOSTK = -1.0
          tagGHOST = 1
          K = ABS(K)
        Endif
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC_Ghost(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          GHOSTI,GHOSTJ,GHOSTK,
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            Call Distance2P_VEC(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          TMPXR, TMPYR, TMPZR)
            ! an external interpolation
            Distance_Real = (TMPXR*TMPX+ TMPYR*TMPY+ TMPZR*TMPZ)
     &                    / Distance
!            If (ABS(Distance_Real) .LT. 1E-5 .OR.
!     &          ABS(Distance/Distance_Real) .GT. 10.0) Then
              PIB2 = P(IB2)
!            Else
!              PIB2 = P(IB) + (P(IB2)-P(IB)) * Distance / Distance_Real
!            Endif
            If (Distance.LE.HSUP .AND. Distance.GT.1E-10
     &         .AND. IB.NE.IB2) Then
                Call KernelGradient(TMPX, TMPY, TMPZ, KX, KY, KZ)
              If (Shepard.EQ.0) Then
                Temp1(IB) = Temp1(IB) -
     &            (P(IB)/Rho(IB)/Rho(IB) + PIB2/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KX
                If (DIMM.EQ.DIMM3D)
     &           Temp2(IB) = Temp2(IB) -
     &            (P(IB)/Rho(IB)/Rho(IB) + PIB2/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KY
                Temp3(IB) = Temp3(IB) -
     &            (P(IB)/Rho(IB)/Rho(IB) + PIB2/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KZ
              Else
                Temp1(IB) = Temp1(IB) -
     &            (P(IB)/Rho(IB)/Rho(IB) + PIB2/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KX
     &            / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     &           Temp2(IB) = Temp2(IB) -
     &            (P(IB)/Rho(IB)/Rho(IB) + PIB2/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KY
     &            / SheCoef(IB)
                Temp3(IB) = Temp3(IB) -
     &            (P(IB)/Rho(IB)/Rho(IB) + PIB2/Rho(IB2)/Rho(IB2))
     &            * Weight(IB2) * Rho(IB2) * KZ
     &            / SheCoef(IB)
              Endif
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
!        if (tagGHOST.GT.0) Then
!          Temp1(IB) = 0.0
!          Temp2(IB) = 0.0
!          Temp3(IB) = 0.0
!        Endif
      Enddo

      End Subroutine

C==================================================================
      Subroutine Calc_VisTerm_Ghost
      include "param.fi"
      include "common.fi"
      INTEGER IB, IJK, II, JJ, KK, DELI, DELJ, DELK, I, J, K, IP, IB2
      REAL TMPX, TMPY, TMPZ, KX, KY, KZ, Distance
      Real GHOSTI, GHOSTJ, GHOSTK, FIB2(3)

      Do IB = 1, NTOT3T
        Temp1(IB) = 0.0
        Temp2(IB) = 0.0
        Temp3(IB) = 0.0
        ! Get Hash Index
        IJK = HASS(IB)
        Call HASHINDEX13(IJK, II, JJ, KK)
        Do DELI = -1, 1
        GHOSTI = 1.0
        I = HASHIND_Ghost(II + DELI)
        If (I .LT. 0) Then
          GHOSTI = -1.0
          I = ABS(I)
        Endif
        Do DELJ = -1, 1
        GHOSTJ = 1.0
        If (DIMM .EQ. DIMM2D .AND. DELJ .NE. 0) Cycle
        J = HASHIND_Ghost(JJ + DELJ)
        If (J .LT. 0) Then
          GHOSTJ = -1.0
          J = ABS(J)
        Endif
        Do DELK = -1, 1
        GHOSTK = 1.0
        K = HASHIND_Ghost(KK + DELK)
        If (K .LT. 0) Then
          GHOSTK = -1.0
          K = ABS(K)
        Endif
          Do IP = 1, HASSNUMS(I,J,K)
            IB2 = HASSIND(IP, I, J, K)
            Call Distance2P_VEC_Ghost(X1(IB2),X2(IB2),X3(IB2),
     &                          X1(IB),X2(IB),X3(IB),
     &                          GHOSTI,GHOSTJ,GHOSTK,
     &                          TMPX, TMPY, TMPZ)
            Distance = SQRT(TMPX*TMPX+ TMPY*TMPY+ TMPZ*TMPZ)
            If (Distance.LE.HSUP) Then
                KX = KernelLaplace(Distance)
                FIB2(1) = U1(IB2)
                FIB2(2) = U2(IB2)
                FIB2(3) = U3(IB2)
                If (GHOSTI.LT.0) FIB2(1) = -FIB2(1) 
                If (GHOSTJ.LT.0) FIB2(2) = -FIB2(2) 
                If (GHOSTK.LT.0) FIB2(3) = -FIB2(3) 
              If (Shepard.EQ.0) Then
                Temp1(IB) = Temp1(IB) +
     &            KX *(FIB2(1)-U1(IB)) * Weight(IB2) * NU
                If (DIMM.EQ.DIMM3D)
     &           Temp2(IB) = Temp2(IB) +
     &            KX *(FIB2(2)-U2(IB)) * Weight(IB2) * NU
                Temp3(IB) = Temp3(IB) +
     &            KX *(FIB2(3)-U3(IB)) * Weight(IB2) * NU
              Else
                Temp1(IB) = Temp1(IB) +
     &            KX *(FIB2(1)-U1(IB)) * Weight(IB2) * NU
     &            / SheCoef(IB)
                If (DIMM.EQ.DIMM3D)
     &           Temp2(IB) = Temp2(IB) +
     &            KX *(FIB2(2)-U2(IB)) * Weight(IB2) * NU
     &            / SheCoef(IB)
                Temp3(IB) = Temp3(IB) +
     &            KX *(FIB2(3)-U3(IB)) * Weight(IB2) * NU
     &            / SheCoef(IB)
              Endif
            Endif
          Enddo
        Enddo
        Enddo
        Enddo
      Enddo

      End Subroutine

